El presente estudio analizará los datos recogidos por una de las estaciones meteorológicas situadas en barcelona con el fin de poder hacer unas predicciones en que condiciones es más probable que precipite sobre Barcelona.
Las predicciones que tienen que ver con el clima están sujetas a gran cantidad de variables que tienen que con la situación geográfica, características físicas del entorno como pueden ser estar cerca del mar o no, estar en un valle o en una zona elevada, etc., pero también con condiciones como puede ser dirección del viento, humedad ambiental, presión, etc.
Aún teniendo en cuenta estos factores, la climatología es impredecible y cambios repentinos en las condiciones atmosféricas, pueden hacer que las predicciones fallen, sobre todo cuando se hacen predicciones de tipo local.
Para hacer el siguiente estudio, tomaremos principalmente los datos de temperatura y si ha habido o no precipitaciones, pero no prescindiremos de otros datos como los meses, dirección viento etc., porque todos estos parámetros pueden hacer variar nuestra predicción final. A priori la hipótesis con la que trabajaremos serà:
* 𝐻0 = la temperatura afecta a las precipitaciones
* H1 = la temperatura no afecta a las precipitaciones.
Los datos los hemos obtenido AEMET (Agencia Española de Meteorología) con la URL, https://opendata.aemet.es/. Los datos poseen la siguiente estrucutura.
* fecha: Fecha del dia (AAAA-MM-DD).
* Indicativo: Indicativo climatológico.
* Nombre: Nombre (ubicación) de la estación.
* provincia: Provincia de la estación.
* altitud: Altitud de la estación en m sobre el nivel del mar.
* tmed: Temperatura media diaria.
* Precipitación diaria de 07 a 07(Ip = inferior a 0,1 mm).
* tmin: Temperatura Mínima del día.
* horatmin: Hora y minuto de la temperatura mínima.
* tmax: Temperatura Máxima del día.
* horatmax: Hora y minuto de la temperatura máxima.
* dir: Dirección de la racha máxima (decenas de grados)
* velmedia: Velocidad media del viento.
* racha: Racha máxima del viento (m/s).
* horaracha: Hora y minuto de la racha máxima.
Hemos descargado los datos desde el 14-06-2017 hasta el 30-06-2021 de la estación meteorológica 0201D de Barcelona . Los datos suministrados por la AEMET anteriores del 14-06-2017, poseen una estructura diferente y menos detallada, por lo que se ha optado por descargar datos posteriores a la fecha indicada hasta el 30 de junio del presente en formato json. El Dataset se ha importado en MySql WorkBench versión 8.0.025 (community).
El presente estudio está organizado de la siguiente manera:
1. Importación de datos y librerías.
2. Exploración y análisis de datos.
2.1 Variables horas y minutos.
2.2 Número de observaciones y valores duplicados.
2.3 La variable 'prec'.
2.4 Variables numéricas.
2.5 Variables Categóricas
2.6 Variables de Temperatura.
3. Outliers/ Imbalanced dataset.
3.1 Imbalanced Dataset.
3.2 Outliers.
4. Análisis no supervisado: clustering.
4.1 Graficar los clusters con PCA.
4.2 Evaluamos los clusters.
4.3 Analizando los cluster.
4.3.1 Analizando las variables categóricas.
4.3.2 Analizando las Variables numéricas.
5. Elegir el modelo.
5.1 Regresión.
5.1.1 Validar el modelo.
5.2 modelos de clasificación.
5.2.1 Validar el rendimiento del modelo
5.2.2 DecisionTreeClassifier
6. Predicción
7. Hipótesis.
!python --version
Python 3.11.0
#Liberias Basicas
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import seaborn as sns
import datetime,time
import re
# importacion librerias SQL
import mysql.connector
#funciones matemáticas
from scipy.stats import ttest_1samp, ttest_ind, f_oneway
#SMOTE
from imblearn.over_sampling import SMOTE
#Clusters
from sklearn.cluster import KMeans
#Evaulación clusters
from sklearn.metrics import silhouette_score
from sklearn.cluster import DBSCAN
#codo
from kneed import KneeLocator
#PCA
from sklearn.decomposition import PCA
#fitter
# from fitter import Fitter, get_common_distributions
#Normalització
from sklearn.preprocessing import LabelEncoder
from sklearn.preprocessing import Normalizer
#Regresión
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
#Classificación
from sklearn.tree import DecisionTreeClassifier
#selection_model
from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import KFold
#metriques regresión
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
from sklearn.metrics import accuracy_score
#Graficar decission tree
import six
import sys
sys.modules['sklearn.externals.six'] = six
from sklearn.tree import export_graphviz
from sklearn.externals.six import StringIO
from IPython.display import Image
# import pydotplus
import warnings
warnings.filterwarnings('ignore')
#Conexión BBDD
mydb = mysql.connector.connect(host="127.0.0.1",user="root",password="Zxcvbn123+-",database="datascience_test")
mydb
#Importacion datos
mycursor = mydb.cursor()
df = pd.read_sql('select * from 20170614_20210630', mydb)
df
Todas las columnas fueron importadas en formato string, para evitar problemas de conversión. Procedemos a dar el formato correcto a aquellas columnas que lo requieran.
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 1455 entries, 0 to 1454 Data columns (total 15 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 fecha 1455 non-null object 1 indicativo 1455 non-null object 2 nombre 1455 non-null object 3 provincia 1455 non-null object 4 altitud 1455 non-null object 5 tmed 1455 non-null object 6 prec 1455 non-null object 7 tmin 1455 non-null object 8 horatmin 1455 non-null object 9 tmax 1455 non-null object 10 horatmax 1455 non-null object 11 dir 1455 non-null object 12 velmedia 1455 non-null object 13 racha 1455 non-null object 14 horaracha 1455 non-null object dtypes: object(15) memory usage: 170.6+ KB
#cambio de object a float
df['tmed']= df['tmed'].str.replace(',','.').astype(float)
df['prec']= df['prec'].str.replace(',','.').astype(float)
df['tmin']= df['tmin'].str.replace(',','.').astype(float)
df['tmax']= df['tmax'].str.replace(',','.').astype(float)
df['dir']= df['dir'].str.replace(',','.').astype(float)
df['racha']= df['racha'].str.replace(',','.').astype(float)
#La velocidad media esta expresada en m/s. La medida temporal que hemos usado es el minuto. Vamos a trasnforar este dato a floa
df['velmedia'] =df['velmedia'].str.replace(',','.').astype(float)
#cambio de object a datetime
df['fecha'] = pd.to_datetime(df['fecha'], format= '%Y/%m/%d')
# Creación de dos columans mes y año.
df['Year'] =df['fecha'].dt.year
df['Month']= df['fecha'].dt.month
df['Day']= df['fecha'].dt.day
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 1455 entries, 0 to 1454 Data columns (total 18 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 fecha 1455 non-null datetime64[ns] 1 indicativo 1455 non-null object 2 nombre 1455 non-null object 3 provincia 1455 non-null object 4 altitud 1455 non-null object 5 tmed 1455 non-null float64 6 prec 1455 non-null float64 7 tmin 1455 non-null float64 8 horatmin 1455 non-null object 9 tmax 1455 non-null float64 10 horatmax 1455 non-null object 11 dir 1455 non-null float64 12 velmedia 1455 non-null float64 13 racha 1455 non-null float64 14 horaracha 1455 non-null object 15 Year 1455 non-null int64 16 Month 1455 non-null int64 17 Day 1455 non-null int64 dtypes: datetime64[ns](1), float64(7), int64(3), object(7) memory usage: 204.7+ KB
Para ayudar a Machine Leanring vamos a convertir las variables horas en minutos.
desglose = df['horatmin'].str.split(':', expand=True)
desglose.columns = ['horatmin_h','horatmin_m']
df= pd.concat([df,desglose], axis=1)
df['horatmin_h'] = pd.to_numeric(df.horatmin_h, errors='coerce')* 60
df['horatmin_m'] = pd.to_numeric(df.horatmin_m, errors='coerce')
df['horatmin_min'] = df['horatmin_h'] + df['horatmin_m']
df =df.drop(['horatmin_h','horatmin_m'], axis=1)
desglose = df['horatmax'].str.split(':', expand=True)
desglose.columns = ['horatmax_h','horatmax_m']
df= pd.concat([df,desglose], axis=1)
df['horatmax_h'] = pd.to_numeric(df.horatmax_h, errors='coerce')* 60
df['horatmax_m'] = pd.to_numeric(df.horatmax_m, errors='coerce')
df['horatmax_min'] = df['horatmax_h'] + df['horatmax_m']
df = df.drop(['horatmax_h','horatmax_m'], axis=1)
desglose = df['horaracha'].str.split(':', expand=True)
desglose.columns = ['horaracha_h','horaracha_m']
df= pd.concat([df,desglose], axis=1)
df['horaracha_h'] = pd.to_numeric(df.horaracha_h, errors='coerce')* 60
df['horaracha_m'] = pd.to_numeric(df.horaracha_m, errors='coerce')
df['horaracha_min'] = df['horaracha_h'] + df['horaracha_m']
df = df.drop(['horaracha_h','horaracha_m'], axis=1)
df.shape
(1455, 21)
#observamos la ausencia de valores nulos
df.isnull().sum()
fecha 0 indicativo 0 nombre 0 provincia 0 altitud 0 tmed 0 prec 0 tmin 0 horatmin 0 tmax 0 horatmax 0 dir 0 velmedia 0 racha 0 horaracha 0 Year 0 Month 0 Day 0 horatmin_min 30 horatmax_min 5 horaracha_min 79 dtype: int64
df[['horatmin_min','horatmax_min','horaracha_min']]=df[['horatmin_min','horatmax_min','horaracha_min']].apply(lambda x: x.fillna(x.mean()),axis=0)
df.describe()
| tmed | prec | tmin | tmax | dir | velmedia | racha | Year | Month | Day | horatmin_min | horatmax_min | horaracha_min | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| count | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 |
| mean | 17.448316 | 1.481512 | 14.364192 | 20.533952 | 24.810997 | 3.534227 | 9.633608 | 2018.978007 | 6.520962 | 15.808247 | 473.593684 | 745.508966 | 757.974564 |
| std | 5.675787 | 6.330925 | 5.935183 | 5.558704 | 21.178533 | 1.518903 | 3.510257 | 1.238231 | 3.434103 | 8.819273 | 401.724058 | 177.900209 | 315.408560 |
| min | 2.700000 | 0.000000 | 0.100000 | 4.600000 | 1.000000 | 0.000000 | 3.900000 | 2017.000000 | 1.000000 | 1.000000 | 10.000000 | 10.000000 | 10.000000 |
| 25% | 12.600000 | 0.000000 | 9.400000 | 15.700000 | 12.000000 | 2.500000 | 7.200000 | 2018.000000 | 4.000000 | 8.000000 | 260.000000 | 650.000000 | 610.000000 |
| 50% | 16.700000 | 0.000000 | 13.600000 | 19.900000 | 21.000000 | 3.300000 | 8.900000 | 2019.000000 | 6.000000 | 16.000000 | 340.000000 | 770.000000 | 800.000000 |
| 75% | 22.450000 | 0.000000 | 19.600000 | 25.450000 | 28.000000 | 4.200000 | 11.400000 | 2020.000000 | 9.000000 | 23.000000 | 440.000000 | 850.000000 | 940.000000 |
| max | 30.000000 | 83.900000 | 27.100000 | 35.200000 | 99.000000 | 18.300000 | 29.400000 | 2021.000000 | 12.000000 | 31.000000 | 1439.000000 | 1439.000000 | 1439.000000 |
El objetivo del estudio es poder hacer análisis sobre las precipitaciones. En consecuencia, a priori, usaremos para las predicciones la variable respuesta 'prec' del dataframe.
df.prec.hist()
<AxesSubplot:>
df.prec
0 0.2
1 0.0
2 0.0
3 0.0
4 0.0
...
1450 0.0
1451 0.0
1452 0.0
1453 0.0
1454 0.0
Name: prec, Length: 1455, dtype: float64
sns.set_theme(style="ticks", palette="pastel")
fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(8, 8))
sns.distplot(
df.prec,
hist = False,
rug = True,
color = "blue",
kde_kws = {'shade': True, 'linewidth': 1},
ax = axes[0]
)
axes[0].set_title("Distribución original", fontsize = 'small')
axes[0].set_xlabel('precipitaciones', fontsize='small')
#axes[0].tick_params(labelsize = 6)
sns.distplot(
np.sqrt(df.prec),
hist = False,
rug = True,
color = "blue",
kde_kws = {'shade': True, 'linewidth': 1},
ax = axes[1]
)
axes[1].set_title("Transformación raíz cuadrada", fontsize = 'small')
axes[1].set_xlabel('sqrt(precipitaciones)', fontsize='small')
axes[1].tick_params(labelsize = 6)
df.skew()
altitud 0.000000 tmed 0.161453 prec 7.304962 tmin 0.147821 tmax 0.128547 dir 2.606666 velmedia 2.381248 racha 1.384965 Year 0.000400 Month -0.002393 Day -0.008205 horatmin_min 1.654524 horatmax_min -0.853173 horaracha_min -0.472763 dtype: float64
distribuciones = ['cauchy', 'chi2', 'expon', 'exponpow', 'gamma',
'norm', 'powerlaw', 'beta', 'logistic']
fitter = Fitter(df.prec, distributions=distribuciones)
fitter.fit()
fitter.summary(Nbest=10, plot=False)
| sumsquare_error | aic | bic | kl_div | |
|---|---|---|---|---|
| expon | 0.353825 | 5745.745876 | -12093.527681 | inf |
| exponpow | 0.653942 | 1205.650080 | -11192.561364 | inf |
| gamma | 0.711120 | 1403.920292 | -11070.599867 | inf |
| logistic | 0.738886 | 6056.698847 | -11022.151867 | inf |
| powerlaw | 0.765632 | 1084.593272 | -10963.132861 | inf |
| beta | 0.772735 | 6863.761516 | -10942.413206 | inf |
| chi2 | 0.775476 | 1339.586925 | -10944.544827 | inf |
| norm | 0.933174 | 6110.042607 | -10682.485264 | inf |
| cauchy | 1.040866 | 10993.477511 | -10523.575373 | inf |
sns.boxplot(x ='Year' ,y ='prec', hue='Year', orient= 'v' ,saturation=0.50, data=df[df.prec > 0],)
<AxesSubplot:xlabel='Year', ylabel='prec'>
La distribución de la variable respuesta es exponencial. Para los modelos lineales generalizados (GML) la distribución ha de ser de la familia de las exponenciales. Es algo a tener en cuenta. En la gràfica de distribuciones podemos observa la presencia de outliers, algo previsible por la naturaleza del dato. Al final del proceso de análisis y exploración procederemos a normalizar estos outliers, porque pueden desvirtuar nuestra predicción respuesta 'prec' del dataframe.
df.select_dtypes(include=['float64', 'int64']).describe()
| tmed | prec | tmin | tmax | dir | velmedia | racha | Year | Month | Day | horatmin_min | horatmax_min | horaracha_min | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| count | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 | 1455.000000 |
| mean | 17.448316 | 1.481512 | 14.364192 | 20.533952 | 24.810997 | 3.534227 | 9.633608 | 2018.978007 | 6.520962 | 15.808247 | 473.593684 | 745.508966 | 757.974564 |
| std | 5.675787 | 6.330925 | 5.935183 | 5.558704 | 21.178533 | 1.518903 | 3.510257 | 1.238231 | 3.434103 | 8.819273 | 401.724058 | 177.900209 | 315.408560 |
| min | 2.700000 | 0.000000 | 0.100000 | 4.600000 | 1.000000 | 0.000000 | 3.900000 | 2017.000000 | 1.000000 | 1.000000 | 10.000000 | 10.000000 | 10.000000 |
| 25% | 12.600000 | 0.000000 | 9.400000 | 15.700000 | 12.000000 | 2.500000 | 7.200000 | 2018.000000 | 4.000000 | 8.000000 | 260.000000 | 650.000000 | 610.000000 |
| 50% | 16.700000 | 0.000000 | 13.600000 | 19.900000 | 21.000000 | 3.300000 | 8.900000 | 2019.000000 | 6.000000 | 16.000000 | 340.000000 | 770.000000 | 800.000000 |
| 75% | 22.450000 | 0.000000 | 19.600000 | 25.450000 | 28.000000 | 4.200000 | 11.400000 | 2020.000000 | 9.000000 | 23.000000 | 440.000000 | 850.000000 | 940.000000 |
| max | 30.000000 | 83.900000 | 27.100000 | 35.200000 | 99.000000 | 18.300000 | 29.400000 | 2021.000000 | 12.000000 | 31.000000 | 1439.000000 | 1439.000000 | 1439.000000 |
sns.set_theme(style="white")
fig, axes = plt.subplots(nrows=6, ncols=2, figsize=(20, 20))
axes = axes.flat
columnas_numeric = df.select_dtypes(include=['float64', 'int64']).columns
columnas_numeric = columnas_numeric.drop('prec')
for i, colum in enumerate(columnas_numeric):
sns.histplot(
data = df,
x = colum,
stat = "count",
kde = True,
color = (list(plt.rcParams['axes.prop_cycle'])*2)[i]["color"],
line_kws= {'linewidth': 2},
alpha = 0.3,
ax = axes[i]
)
axes[i].set_title(colum, fontsize = 25, fontweight = "bold")
axes[i].tick_params(labelsize = 20)
axes[i].set_xlabel("")
for i in [11]:
fig.delaxes(axes[i])
fig.tight_layout()
plt.subplots_adjust(top = 0.9)
fig.suptitle('Distribución variables numéricas', fontsize = 25, fontweight = "bold");
relación variables numéricas con variable respuesta.
dfcorr = df.corr()
pd.set_option('float_format', '{:.2f}'.format)
f = plt.subplots(figsize=(11, 9))
sns.heatmap(dfcorr,annot=True, annot_kws={"size":11}, square=True, linewidths=.5, cbar_kws={"shrink": .5}, robust=False, center=0, vmin=-1, vmax=1)
plt.show()
dfvar = df.var()
dfvar
tmed 32.21 prec 40.08 tmin 35.23 tmax 30.90 dir 448.53 velmedia 2.31 racha 12.32 Year 1.53 Month 11.79 Day 77.78 horatmin_min 161382.22 horatmax_min 31648.48 horaracha_min 99482.56 dtype: float64
La correlación entre las variables numéricas en relación a prec no son muy buenas. Difícilmente podremos hacer una modelo de regresión con buenos resultados.
dfcorr.sort_values(by='prec')
| tmed | prec | tmin | tmax | dir | velmedia | racha | Year | Month | Day | horatmin_min | horatmax_min | horaracha_min | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| horatmax_min | -0.17 | -0.12 | -0.18 | -0.14 | 0.03 | 0.01 | 0.04 | 0.01 | -0.04 | 0.00 | -0.18 | 1.00 | 0.07 |
| dir | -0.06 | -0.09 | -0.06 | -0.04 | 1.00 | -0.15 | -0.13 | -0.02 | 0.04 | -0.02 | -0.05 | 0.03 | -0.01 |
| tmax | 0.99 | -0.05 | 0.95 | 1.00 | -0.04 | 0.03 | -0.14 | -0.20 | 0.30 | 0.03 | -0.13 | -0.14 | 0.11 |
| tmed | 1.00 | -0.03 | 0.99 | 0.99 | -0.06 | 0.07 | -0.13 | -0.19 | 0.30 | 0.03 | -0.10 | -0.17 | 0.11 |
| Year | -0.19 | -0.02 | -0.18 | -0.20 | -0.02 | 0.06 | 0.03 | 1.00 | -0.35 | -0.02 | 0.00 | 0.01 | 0.03 |
| tmin | 0.99 | -0.01 | 1.00 | 0.95 | -0.06 | 0.11 | -0.11 | -0.18 | 0.29 | 0.03 | -0.08 | -0.18 | 0.10 |
| Day | 0.03 | 0.00 | 0.03 | 0.03 | -0.02 | -0.01 | -0.04 | -0.02 | 0.01 | 1.00 | -0.02 | 0.00 | -0.01 |
| Month | 0.30 | 0.03 | 0.29 | 0.30 | 0.04 | -0.12 | -0.07 | -0.35 | 1.00 | 0.01 | 0.03 | -0.04 | 0.02 |
| horaracha_min | 0.11 | 0.06 | 0.10 | 0.11 | -0.01 | 0.13 | 0.05 | 0.03 | 0.02 | -0.01 | -0.11 | 0.07 | 1.00 |
| horatmin_min | -0.10 | 0.15 | -0.08 | -0.13 | -0.05 | 0.08 | 0.24 | 0.00 | 0.03 | -0.02 | 1.00 | -0.18 | -0.11 |
| velmedia | 0.07 | 0.21 | 0.11 | 0.03 | -0.15 | 1.00 | 0.64 | 0.06 | -0.12 | -0.01 | 0.08 | 0.01 | 0.13 |
| racha | -0.13 | 0.25 | -0.11 | -0.14 | -0.13 | 0.64 | 1.00 | 0.03 | -0.07 | -0.04 | 0.24 | 0.04 | 0.05 |
| prec | -0.03 | 1.00 | -0.01 | -0.05 | -0.09 | 0.21 | 0.25 | -0.02 | 0.03 | 0.00 | 0.15 | -0.12 | 0.06 |
Las Variables categóricas que había en el dataset o no tienen relevancia porque, o son constantes o bien las hemos transformado a numericas y analizado en apartados anteriores.
Visualizamos la temperatura media, max y minima porque parecen tener relación con las precipitaciones
plt.style.use('ggplot')
fig = plt.figure(figsize=(10, 10), dpi=80)
ax = fig.add_subplot(3,1,1)
x = df['tmed'].value_counts().keys()
y = df['tmed'].value_counts().values
plt.title('Temperaturas medias Barcelona')
plt.bar(x, y)
ax1 = fig.add_subplot(3,1,2) #3,1,2
x = df['tmin'].value_counts().keys()
y = df['tmin'].value_counts().values
plt.bar(x, y)
plt.title('Temperaturas mínimas Barcelona')
ax2 = fig.add_subplot(3,1,3) # 3,1,3
x = df['tmax'].value_counts().keys()
y = df['tmax'].value_counts().values
plt.bar(x, y)
plt.title('Temperaturas máximas Barcelona')
fig.tight_layout()
plt.show
<function matplotlib.pyplot.show(close=None, block=None)>
Por la naturaleza del dato es previsible, que el número de días que no ha llovido sea superior al número de días que no lo ha hecho. Este hecho produce que el dataset esté descompensado. (Imbalanced Dataset). Aplicaremos el algoritmo SMOTE para generar datos “sintéticos”, que nos igualen el número de días que ha llovido, respecto a los que no lo ha hecho. De esta manera nuestras predicciones serán más ajustadas.
df_imb = pd.DataFrame(df)
df_imb['prec2'] = df_imb['prec'].apply(lambda x : 0 if x == 0 else 1).astype('int64')
fig = plt.figure(figsize=(4,4))
ax = sns.factorplot('prec2',data=df_imb,kind='count',aspect=2)
plt.show()
<Figure size 288x288 with 0 Axes>
df_imb = df_imb[['prec2','tmed', 'tmin', 'tmax', 'dir', 'velmedia', 'racha',
'horatmin_min', 'horatmax_min', 'horaracha_min']]
df_imb
| prec2 | tmed | tmin | tmax | dir | velmedia | racha | horatmin_min | horatmax_min | horaracha_min | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 25.00 | 22.00 | 28.00 | 20.00 | 4.20 | 8.90 | 280.00 | 860.00 | 870.00 |
| 1 | 0 | 24.80 | 22.50 | 27.20 | 20.00 | 2.80 | 6.90 | 50.00 | 930.00 | 840.00 |
| 2 | 0 | 26.30 | 23.30 | 29.30 | 24.00 | 4.40 | 9.20 | 200.00 | 1020.00 | 1000.00 |
| 3 | 0 | 26.20 | 24.30 | 28.00 | 7.00 | 3.10 | 5.60 | 310.00 | 690.00 | 370.00 |
| 4 | 0 | 25.20 | 22.60 | 27.80 | 20.00 | 4.20 | 6.90 | 1439.00 | 710.00 | 830.00 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1450 | 0 | 22.00 | 19.60 | 24.30 | 20.00 | 5.60 | 12.50 | 240.00 | 470.00 | 870.00 |
| 1451 | 0 | 22.80 | 20.30 | 25.20 | 20.00 | 5.00 | 10.00 | 170.00 | 940.00 | 740.00 |
| 1452 | 0 | 23.40 | 20.80 | 25.90 | 20.00 | 2.80 | 6.70 | 300.00 | 840.00 | 650.00 |
| 1453 | 0 | 23.40 | 21.50 | 25.20 | 14.00 | 4.20 | 11.10 | 230.00 | 560.00 | 1420.00 |
| 1454 | 0 | 23.40 | 21.60 | 25.20 | 4.00 | 4.70 | 9.70 | 250.00 | 660.00 | 1070.00 |
1455 rows × 10 columns
Xsmt = df_imb.loc[:, df_imb.columns != 'prec2'].values
ysmt = df_imb.prec2
sm = SMOTE(sampling_strategy='auto', k_neighbors=1, random_state=100)
Xbal, ybla = sm.fit_resample(Xsmt, ysmt)
dfsmote = pd.concat([pd.DataFrame(Xbal),pd.DataFrame(ybla)], axis = 1)
dfsmote.columns= ['tmed', 'tmin', 'tmax', 'dir', 'velmedia', 'racha',
'horatmin_min', 'horatmax_min', 'horaracha_min','prec2']
dfsmote
| tmed | tmin | tmax | dir | velmedia | racha | horatmin_min | horatmax_min | horaracha_min | prec2 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 25.00 | 22.00 | 28.00 | 20.00 | 4.20 | 8.90 | 280.00 | 860.00 | 870.00 | 1 |
| 1 | 24.80 | 22.50 | 27.20 | 20.00 | 2.80 | 6.90 | 50.00 | 930.00 | 840.00 | 0 |
| 2 | 26.30 | 23.30 | 29.30 | 24.00 | 4.40 | 9.20 | 200.00 | 1020.00 | 1000.00 | 0 |
| 3 | 26.20 | 24.30 | 28.00 | 7.00 | 3.10 | 5.60 | 310.00 | 690.00 | 370.00 | 0 |
| 4 | 25.20 | 22.60 | 27.80 | 20.00 | 4.20 | 6.90 | 1439.00 | 710.00 | 830.00 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2265 | 17.12 | 15.03 | 19.13 | 8.78 | 2.83 | 8.89 | 55.56 | 785.53 | 942.22 | 1 |
| 2266 | 20.90 | 16.94 | 24.94 | 29.91 | 3.12 | 10.32 | 594.29 | 829.05 | 761.97 | 1 |
| 2267 | 22.78 | 20.34 | 25.14 | 19.14 | 3.14 | 8.45 | 1110.00 | 840.00 | 743.74 | 1 |
| 2268 | 13.58 | 12.05 | 15.19 | 13.58 | 7.76 | 15.42 | 459.94 | 755.06 | 453.56 | 1 |
| 2269 | 16.80 | 14.52 | 19.08 | 7.23 | 3.76 | 9.54 | 1402.21 | 655.47 | 48.42 | 1 |
2270 rows × 10 columns
fig = plt.figure(figsize=(4,4))
ax = sns.factorplot('prec2',data=dfsmote,kind='count',aspect=2)
plt.show()
<Figure size 288x288 with 0 Axes>
Ahora ya tenemos el dataset compensado, hemos pasado de tener 1455 registros iniciales a 2270.
Como hemos visto en el apartado 2.4 sobre las variables numéricas, tenemos puntas de temperatura, velocidad de vientos que son excepcionales y que pueden alterar la predicciones, por ello vamos a procesar los datos con el fin de eliminar estas medidas que estan fuera de lo que seria una distribución normal (outliers).
df2 = pd.DataFrame(dfsmote)
Q1 = df2.quantile(0.10)
Q3 = df2.quantile(0.90)
IQR = Q3 - Q1
df2.shape
(2270, 10)
df2 = df2[~((df2 < (Q1 - 1.5 * IQR)) |(df2 > (Q3 + 1.5 * IQR))).any(axis=1)]
df2.shape, df.shape
((2128, 10), (1455, 22))
Hemos pasado de tener 2270 registros a 2128. El Dataframe df2, lo reservaremos para aplicar los modelos predictivos más adelante.
Crearemos un dataframe con las variables normalizadas, seleccionando las variables que necesitemos para el clustering. Omitiremos aquellas columnas que no aportan información como son indicativo, nombre, provincia, altitud que son valores constantes.
scalerN = Normalizer()
#normalizamos los datos
df_s = scalerN.fit_transform(df2)
df_s
array([[1.99060054e-02, 1.75172847e-02, 2.22947260e-02, ...,
6.84766584e-01, 6.92728987e-01, 7.96240214e-04],
[1.97591273e-02, 1.79266276e-02, 2.16713009e-02, ...,
7.40967273e-01, 6.69260763e-01, 0.00000000e+00],
[1.82218754e-02, 1.61433345e-02, 2.03004163e-02, ...,
7.06703912e-01, 6.92846973e-01, 0.00000000e+00],
...,
[1.44304194e-02, 1.28796648e-02, 1.59232788e-02, ...,
5.32021970e-01, 4.71057513e-01, 6.33359488e-04],
[1.36610341e-02, 1.21220822e-02, 1.52789735e-02, ...,
7.59471153e-01, 4.56207411e-01, 1.00583990e-03],
[1.08461500e-02, 9.37620969e-03, 1.23160903e-02, ...,
4.23174827e-01, 3.12616126e-02, 6.45604167e-04]])
sse = []
for k in range(1, 10):
kmeans = KMeans(n_clusters=k)
kmeans.fit(df_s)
sse.append(kmeans.inertia_)
plt.style.use("fivethirtyeight")
plt.plot(range(1, 10), sse)
plt.xticks(range(1, 10))
plt.xlabel("Numero de clusters")
plt.ylabel("SSE")
plt.show()
kl = KneeLocator(range(1, 10), sse,
curve="convex",
direction="decreasing")
print('El numero de clusters debe ser :', kl.elbow)
El numero de clusters debe ser : 3
clustering = KMeans(n_clusters=kl.elbow,max_iter=300)
clustering.fit(df_s)
KMeans(n_clusters=3)
# Añadimos el número de cluster al dataframe original para su analisis.
df2['Kmeans_Clustering']= clustering.labels_
pca = PCA(n_components=2)
pca1 = pca.fit_transform (df_s)
df_pca = pd.DataFrame(pca1, columns=['componente1', 'componente2'])
df_pca['Kmeans_Clustering']= clustering.labels_
df_pca
| componente1 | componente2 | Kmeans_Clustering | |
|---|---|---|---|
| 0 | -0.22 | 0.02 | 0 |
| 1 | -0.37 | 0.10 | 0 |
| 2 | -0.29 | 0.05 | 0 |
| 3 | -0.02 | 0.26 | 2 |
| 4 | 0.45 | -0.12 | 1 |
| ... | ... | ... | ... |
| 2123 | -0.38 | -0.04 | 0 |
| 2124 | 0.03 | 0.02 | 0 |
| 2125 | 0.33 | -0.01 | 1 |
| 2126 | 0.07 | 0.20 | 2 |
| 2127 | 0.75 | 0.18 | 1 |
2128 rows × 3 columns
fig = plt.figure(figsize=(4,5))
color_t = np.array(['green','blue','yellow'])
ax = fig.add_subplot(1,1,1)
ax.set_xlabel('Principal Component 1', fontsize = 10)
ax.set_ylabel('Principal Component 2', fontsize = 10)
ax.set_title('2 component PCA', fontsize = 12)
ax.scatter(x = df_pca.componente1, y= df_pca.componente2,
c = color_t[df_pca.Kmeans_Clustering], s=4)
<matplotlib.collections.PathCollection at 0x150102bc5e0>
Observamos que tenemos tres clusters bien diferenciados.
print('Evaulación Kmeans_silhouette:', silhouette_score(df_s, clustering.labels_).round(1))
Evaulación Kmeans_silhouette: 0.6
Silhouette se refiere a un método de interpretación y validación de la coherencia dentro del análisis de grupos, que oscila entre -1 y 1, por lo que podemos afirmar que la configuración es correcta.
Comenzaremos observando algunos valores de los diferentes clusters. En primer lugar una vista general de la información registrada en cada cluster.
df2.info()
<class 'pandas.core.frame.DataFrame'> Int64Index: 2128 entries, 0 to 2269 Data columns (total 11 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 tmed 2128 non-null float64 1 tmin 2128 non-null float64 2 tmax 2128 non-null float64 3 dir 2128 non-null float64 4 velmedia 2128 non-null float64 5 racha 2128 non-null float64 6 horatmin_min 2128 non-null float64 7 horatmax_min 2128 non-null float64 8 horaracha_min 2128 non-null float64 9 prec2 2128 non-null int64 10 Kmeans_Clustering 2128 non-null int32 dtypes: float64(9), int32(1), int64(1) memory usage: 191.2 KB
df2['Kmeans_Clustering'] = df2.Kmeans_Clustering.astype('int64')
print('Componentes por cluster :', df2.groupby('Kmeans_Clustering').size())
Componentes por cluster : Kmeans_Clustering 0 1255 1 514 2 359 dtype: int64
print('Mediana Precipitaciones por Cluster: ', df2.groupby('Kmeans_Clustering')['prec2'].median())
print('Media Precipitaciones por Cluster:', df2.groupby('Kmeans_Clustering')['prec2'].mean())
print('Media Precipitaciones por Cluster:', df2.groupby('Kmeans_Clustering')['prec2'].max())
Mediana Precipitaciones por Cluster: Kmeans_Clustering 0 0 1 1 2 0 Name: prec2, dtype: int64 Media Precipitaciones por Cluster: Kmeans_Clustering 0 0.43 1 0.77 2 0.41 Name: prec2, dtype: float64 Media Precipitaciones por Cluster: Kmeans_Clustering 0 1 1 1 2 1 Name: prec2, dtype: int64
df2.groupby(['Kmeans_Clustering'])['prec2'].sum().round(2)
Kmeans_Clustering 0 534 1 395 2 147 Name: prec2, dtype: int64
No hay variables categoricas en el Dataframe seleccionado.
Analizaremos cada cluster por separado.
df_C =pd.DataFrame(df2.select_dtypes(include=['float','int64']))
df_C_estadistica = pd.DataFrame()
#df_C_estadistica['Columna']= df_C.columns
df_C0= df2[df_C.Kmeans_Clustering==0]
df_C1= df2[df_C.Kmeans_Clustering==1]
df_C2= df2[df_C.Kmeans_Clustering==2]
df_C_estadistica['mean cluster 0'] = df_C0.mean()
df_C_estadistica['std cluster 0'] = df_C0.std()
df_C_estadistica['mean cluster 1'] = df_C1.mean()
df_C_estadistica['std cluster 1'] = df_C1.std()
df_C_estadistica['mean cluster 2'] = df_C2.mean()
df_C_estadistica['std cluster 2'] = df_C2.std()
df_C_estadistica.round(2)
| mean cluster 0 | std cluster 0 | mean cluster 1 | std cluster 1 | mean cluster 2 | std cluster 2 | |
|---|---|---|---|---|---|---|
| tmed | 17.62 | 5.42 | 16.95 | 5.08 | 15.95 | 5.33 |
| tmin | 14.52 | 5.64 | 14.18 | 4.98 | 12.89 | 5.80 |
| tmax | 20.72 | 5.33 | 19.71 | 5.30 | 19.00 | 5.01 |
| dir | 19.60 | 8.07 | 17.22 | 9.07 | 20.87 | 11.84 |
| velmedia | 3.58 | 1.21 | 4.04 | 1.54 | 3.41 | 1.64 |
| racha | 9.69 | 2.92 | 12.74 | 3.85 | 9.41 | 3.39 |
| horatmin_min | 301.28 | 109.21 | 1278.02 | 240.87 | 325.64 | 137.52 |
| horatmax_min | 767.42 | 154.19 | 602.39 | 266.47 | 805.13 | 141.66 |
| horaracha_min | 924.77 | 199.83 | 692.49 | 358.72 | 308.85 | 186.45 |
| prec2 | 0.43 | 0.49 | 0.77 | 0.42 | 0.41 | 0.49 |
| Kmeans_Clustering | 0.00 | 0.00 | 1.00 | 0.00 | 2.00 | 0.00 |
Observamos que tanto la media como la desviación estándar de las observaciones son bastante similares. En cuanto a la media de las precipitaciones en el cluster 2 es de 4.01 mientras que en los clusters 0 y 1 son 0,43 y 0,77 respectivamente.
Por una lado, en el apartado 2.3 variable respuesta, observamos que la distribución de la variable respuesta era exponencial. Este tipo de distribución se adapta como norma general a los modelos de regresión, pero no había buena correlación entre nuestra variable respuesta “prec2” y con el resto de variables. Es previsible que una regresión lineal con los datos que tenemos, no sea una buena elección.
df2
| tmed | tmin | tmax | dir | velmedia | racha | horatmin_min | horatmax_min | horaracha_min | prec2 | Kmeans_Clustering | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 25.00 | 22.00 | 28.00 | 20.00 | 4.20 | 8.90 | 280.00 | 860.00 | 870.00 | 1 | 0 |
| 1 | 24.80 | 22.50 | 27.20 | 20.00 | 2.80 | 6.90 | 50.00 | 930.00 | 840.00 | 0 | 0 |
| 2 | 26.30 | 23.30 | 29.30 | 24.00 | 4.40 | 9.20 | 200.00 | 1020.00 | 1000.00 | 0 | 0 |
| 3 | 26.20 | 24.30 | 28.00 | 7.00 | 3.10 | 5.60 | 310.00 | 690.00 | 370.00 | 0 | 2 |
| 4 | 25.20 | 22.60 | 27.80 | 20.00 | 4.20 | 6.90 | 1439.00 | 710.00 | 830.00 | 0 | 1 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2265 | 17.12 | 15.03 | 19.13 | 8.78 | 2.83 | 8.89 | 55.56 | 785.53 | 942.22 | 1 | 0 |
| 2266 | 20.90 | 16.94 | 24.94 | 29.91 | 3.12 | 10.32 | 594.29 | 829.05 | 761.97 | 1 | 0 |
| 2267 | 22.78 | 20.34 | 25.14 | 19.14 | 3.14 | 8.45 | 1110.00 | 840.00 | 743.74 | 1 | 1 |
| 2268 | 13.58 | 12.05 | 15.19 | 13.58 | 7.76 | 15.42 | 459.94 | 755.06 | 453.56 | 1 | 2 |
| 2269 | 16.80 | 14.52 | 19.08 | 7.23 | 3.76 | 9.54 | 1402.21 | 655.47 | 48.42 | 1 | 1 |
2128 rows × 11 columns
#df2 = df2[['prec','tmed','tmin', 'tmax']]
df2 = df2[['prec2','tmed', 'tmin', 'tmax', 'dir', 'velmedia', 'racha',
'horatmin_min', 'horatmax_min', 'horaracha_min']]
yr = np.array(df2.prec2).reshape(-1,1)
Xr = np.array(df2.drop(['prec2'],axis=1))
X_train, X_test, y_train, y_test = train_test_split(Xr,yr, test_size=0.3,random_state=4 )
X_train= scalerN.fit_transform(X_train)
#recuperamos las variable X y y generadas en el apartado 3.
models=[LinearRegression(), Ridge(), Lasso(),DecisionTreeRegressor(),RandomForestRegressor(random_state=32)]
for model in models:
model.fit(X_train,y_train)
prediccions = model.predict(X_test)
print(type(model).__name__)
print(" MAE", mean_absolute_error(y_test,prediccions))
print(" RMSE", (mean_squared_error(y_test,prediccions)))
print(" R2", r2_score(y_test,prediccions))
LinearRegression MAE 244.1340755145745 RMSE 101683.6301705869 R2 -406854.0862801001 Ridge MAE 376.6891279911446 RMSE 161559.8815795461 R2 -646430.0867853195 Lasso MAE 0.4999248532009909 RMSE 0.249943909445553 R2 -7.19945979192893e-05 DecisionTreeRegressor MAE 0.49139280125195617 RMSE 0.49139280125195617 R2 -0.9661538461538464 RandomForestRegressor MAE 0.49871674491392803 RMSE 0.2523067292644757 R2 -0.00952607545320916
Hemos testado varios modelos de regressión y en ninguno de ellos el resultado es bueno.
df2
| prec2 | tmed | tmin | tmax | dir | velmedia | racha | horatmin_min | horatmax_min | horaracha_min | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 25.00 | 22.00 | 28.00 | 20.00 | 4.20 | 8.90 | 280.00 | 860.00 | 870.00 |
| 1 | 0 | 24.80 | 22.50 | 27.20 | 20.00 | 2.80 | 6.90 | 50.00 | 930.00 | 840.00 |
| 2 | 0 | 26.30 | 23.30 | 29.30 | 24.00 | 4.40 | 9.20 | 200.00 | 1020.00 | 1000.00 |
| 3 | 0 | 26.20 | 24.30 | 28.00 | 7.00 | 3.10 | 5.60 | 310.00 | 690.00 | 370.00 |
| 4 | 0 | 25.20 | 22.60 | 27.80 | 20.00 | 4.20 | 6.90 | 1439.00 | 710.00 | 830.00 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 2265 | 1 | 17.12 | 15.03 | 19.13 | 8.78 | 2.83 | 8.89 | 55.56 | 785.53 | 942.22 |
| 2266 | 1 | 20.90 | 16.94 | 24.94 | 29.91 | 3.12 | 10.32 | 594.29 | 829.05 | 761.97 |
| 2267 | 1 | 22.78 | 20.34 | 25.14 | 19.14 | 3.14 | 8.45 | 1110.00 | 840.00 | 743.74 |
| 2268 | 1 | 13.58 | 12.05 | 15.19 | 13.58 | 7.76 | 15.42 | 459.94 | 755.06 | 453.56 |
| 2269 | 1 | 16.80 | 14.52 | 19.08 | 7.23 | 3.76 | 9.54 | 1402.21 | 655.47 | 48.42 |
2128 rows × 10 columns
df2 = df2[['prec2','tmed', 'tmin', 'tmax', 'dir', 'velmedia', 'racha',
'horatmin_min', 'horatmax_min', 'horaracha_min']]
#feature_cols = df3.columns
feature_cols = ['tmed', 'tmin', 'tmax', 'dir', 'velmedia', 'racha',
'horatmin_min', 'horatmax_min', 'horaracha_min']
yc = np.array(df2.prec2)
Xc = np.array(df2.drop(['prec2'],axis=1))
X_train, X_test, y_train, y_test = train_test_split(Xc, yc, test_size=0.3, random_state=1)
#X_train= scalerN.fit_transform(X_train)
# Create Decision Tree classifer object
clf = DecisionTreeClassifier()
# Train Decision Tree Classifer
clf = clf.fit(X_train,y_train)
#Predict the response for test dataset
y_pred = clf.predict(X_test)
kf = KFold(n_splits=5)
score =clf.score(X_train,y_train)
scores = cross_val_score(clf, X_train, y_train, cv=kf, scoring="accuracy")
score_pred = accuracy_score(y_test, y_pred)
print("Metricas cross_validation", scores)
print("Metrica del modelo", score)
print("Media de cross_validation", scores.mean())
print("Metrica en Test", score_pred)
Metricas cross_validation [0.77181208 0.74832215 0.77852349 0.77516779 0.74074074] Metrica del modelo 1.0 Media de cross_validation 0.7629132488192891 Metrica en Test 0.7996870109546166
La comprobación del rendimiento es buena por lo que proseguimos con nuestro modelo.
dot_data = StringIO()
export_graphviz(clf, out_file=dot_data,
filled=True, rounded=True,
special_characters=True,feature_names = feature_cols,class_names=['0','1'])
export_graphviz
<function sklearn.tree._export.export_graphviz(decision_tree, out_file=None, *, max_depth=None, feature_names=None, class_names=None, label='all', filled=False, leaves_parallel=False, impurity=True, node_ids=False, proportion=False, rotate=False, rounded=False, special_characters=False, precision=3)>
graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
graph.write_png('predicción_bcn2.png')
True
Image(graph.create_png())
acc_decision_tree = round(clf.score(X_train, y_train) * 100, 2)
print(acc_decision_tree)
100.0
df_lr=pd.DataFrame({'Actual':y_test, 'Predicted':y_pred}, index=list(range(0, 639)))
df_lr.round(3)
| Actual | Predicted | |
|---|---|---|
| 0 | 1 | 0 |
| 1 | 1 | 1 |
| 2 | 1 | 0 |
| 3 | 0 | 0 |
| 4 | 0 | 0 |
| ... | ... | ... |
| 634 | 1 | 1 |
| 635 | 0 | 1 |
| 636 | 0 | 0 |
| 637 | 0 | 1 |
| 638 | 1 | 1 |
639 rows × 2 columns
Vamos a analizar el resultado de nuestra hipótesis.
maximo = df2['tmax'][df2['prec2'] == 1]
minimo = df2['tmin'][df2['prec2'] == 1]
medio = df2['tmed'][df2['prec2'] == 1]
maximo,minimo,medio
(0 28.00
11 27.90
14 30.50
15 27.00
39 27.90
...
2265 19.13
2266 24.94
2267 25.14
2268 15.19
2269 19.08
Name: tmax, Length: 1076, dtype: float64,
0 22.00
11 21.20
14 20.20
15 19.50
39 22.00
...
2265 15.03
2266 16.94
2267 20.34
2268 12.05
2269 14.52
Name: tmin, Length: 1076, dtype: float64,
0 25.00
11 24.60
14 25.40
15 23.20
39 25.00
...
2265 17.12
2266 20.90
2267 22.78
2268 13.58
2269 16.80
Name: tmed, Length: 1076, dtype: float64)
# test ANOVA
alfa = 0.05
stat, p = f_oneway(maximo, minimo, medio)
print(f'F-statistic = {stat:.3f}\np-value = {p:.3f}')
print('No rechazamos H0') if p > alfa else print('Rechazamos H0')
F-statistic = 347.194 p-value = 0.000 Rechazamos H0
Aplicamos el método Nova para comprobar nuestra hipótesis, que queda rechazada.